home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_300 / 315_01 / fft.c < prev    next >
C/C++ Source or Header  |  1990-05-14  |  14KB  |  566 lines
  1.  
  2. /* int_scale(), dbl_scale(), and long_scale() moved from herc.c */
  3. /* 10/89, because herc.c is being discontinued, and msherc.com is */
  4. /* being used for hercules graphcis support.  T Clune */
  5.  
  6. /* fft.c is a file of functions to support fourier transforms. */
  7. /* Written 3/89 by T Clune.  Copyright (c) 1989, Eye Research Institute, */
  8. /* 20 Staniford St, Boston, MA 02114.  All Rights Reserved. */
  9.  
  10. #include <stdio.h>
  11. #include <stdlib.h>
  12. #include <math.h>
  13. #include <malloc.h>
  14. #include <conio.h>
  15. #include "fft.h"
  16.  
  17.  
  18.  
  19. /* data_multiply() multiplies two data files of equal length together */
  20.  
  21. void data_multiply(file1, file2)
  22. char file1[], file2[];
  23. {
  24.     FILE *f1, *f2, *f3;
  25.     char string[80];
  26.     int i, numpts1, numpts2;
  27.     double d1, d2, *d_out, minval, maxval;
  28.  
  29.  
  30.     f1=fopen(file1,"r");
  31.     if(f1==NULL)
  32.     {
  33.     printf("Error opening data file.  Exiting routine.\n");
  34.     printf("Press any key to continue...\n");
  35.     mgetch();
  36.     return;
  37.     }
  38.  
  39.     f2=fopen(file2,"r");
  40.     if(f2==NULL)
  41.     {
  42.     printf("Error opening data file.  Exiting routine.\n");
  43.     printf("Press any key to continue...\n");
  44.     mgetch();
  45.     return;
  46.     }
  47.  
  48.     fscanf(f1,"%d %*lf %*lf", &numpts1);
  49.     fscanf(f2,"%d %*lf %*lf", &numpts2);
  50.     if(numpts1 != numpts2)
  51.     {
  52.     printf("Different size input files.  Exiting routine.\n");
  53.     printf("Press any key to continue...\n");
  54.     mgetch();
  55.     return;
  56.     }
  57.  
  58.     while(kbhit()) ;
  59.     printf("Enter filespec for output file\n");
  60.     mgets(string);
  61.     f3=fopen(string,"w");
  62.  
  63.     if(f3==NULL)
  64.     {
  65.     printf("Error opening output file.  Exiting routine.\n");
  66.     printf("Press any key to continue...\n");
  67.     mgetch();
  68.     return;
  69.     }
  70.     d_out=(double *)malloc(sizeof(double)*numpts1);
  71.     if(d_out==NULL)
  72.     {
  73.     printf("malloc() error in data_multiply(). Program aborting.\n");
  74.     exit(-1);
  75.     }
  76.  
  77.  
  78.     for(i=0;i<numpts1;i++)
  79.     {
  80.     fscanf(f1,"%lf", &d1);
  81.     fscanf(f2,"%lf",&d2);
  82.     d_out[i]=d1*d2;
  83.     }
  84.  
  85.     minmax(d_out,numpts1,&minval,&maxval);
  86.     fprintf(f3,"%d %f %f\n", numpts1, minval, maxval);
  87.  
  88.     for(i=0;i<numpts1;i++)
  89.     fprintf(f3,"%f\n",d_out[i]);
  90.  
  91.     fclose(f1);
  92.     fclose(f2);
  93.     fclose(f3);
  94.     free(d_out);
  95.  
  96. }
  97.  
  98.  
  99.  
  100.  
  101. /* dbl_scale() scales an extended precision floating point number */
  102. /* to an integer value.  It is useful for scaling floating point  */
  103. /* numbers for graphing on the screen.                            */
  104. /* moved from herc.c 10/89 by T Clune */
  105.  
  106.  
  107. int dbl_scale(val, min, max, out_min, out_max)
  108.  
  109. double val, min, max;  /* val= number to be scaled; min= smallest possible */
  110.     /* value for val; max=largest possible value */
  111. int out_min, out_max;  /* out_min=integer to return if val=min; */
  112.         /* out_max=integer to return if val=max.  NOTE THAT out_min */
  113.         /* MAY BE LARGER THAN out_max.  If you want to change screen */
  114.         /* values from 4th quadrant to 1st quadrant, the y-value */
  115.         /* out_min would be larger than the y-value out_max, while */
  116.         /* x-value out_max would be larger than the x-value out_min */
  117.  
  118. {
  119.  
  120.     int out_val, largest, smallest; /* val returned, larger of out_min */
  121.         /* and out_max, smaller of the same */
  122.  
  123.     double scale_factor;    /* ratio of in and out ranges */
  124.         /* scale the input value */
  125.  
  126.     scale_factor = (out_max - out_min) / (max - min);
  127.     out_val = (int)((val - min) * scale_factor + out_min);
  128.  
  129.  
  130.     if(out_min>out_max)     /* check for max/min inversion */
  131.     {
  132.     smallest = out_max;
  133.     largest = out_min;
  134.  
  135.     }
  136.     else
  137.     {
  138.     smallest = out_min;
  139.     largest = out_max;
  140.     }
  141.  
  142.     if(out_val > largest)      /* check for out-of-bounds values */
  143.     out_val = largest;     /* since val is not checked to insure */
  144.     if(out_val < smallest)      /* that it is between min and max, this */
  145.     out_val = smallest;     /* is necessary.  By checking at the last */
  146.                 /* step, any arithmetic round-off accumulation */
  147.                 /* is also trapped */
  148.  
  149.     return out_val;          /* return the scaled value */
  150.  
  151. }
  152.  
  153.  
  154.  
  155.  
  156.  
  157. /* fft() accepts the time-domain data x,y (n points, where n is a power of 2) */
  158. /* and converts to frequency domain rectangular data if flag==FORWARD, */
  159. /* or accepts frequency-domain data in x,y and converts to time domain if */
  160. /* flag==INVERSE.  FORWARD and INVERSE are defined in fft.h */
  161.  
  162. void fft(x, y, n, flag)
  163. double x[], y[];
  164. int n, flag;
  165. {
  166.    int maxpower, arg, pnt0, pnt1;
  167.    int i, j=0, a, b, k;
  168.    double prodreal, prodimag, harm, temp;
  169.    double *cos_coef;
  170.    double *sin_coef;
  171.  
  172.    cos_coef=(double*)calloc(n,sizeof(double));
  173.    sin_coef=(double*)calloc(n,sizeof(double));
  174.  
  175.    if((cos_coef==NULL) || (sin_coef==NULL))
  176.    {
  177.     printf("malloc() error in fft function\n");
  178.     if(cos_coef != NULL)
  179.         free(cos_coef);
  180.     if(sin_coef != NULL)
  181.         free(sin_coef);
  182.     return;
  183.     }
  184.  
  185.    if(flag==INVERSE)
  186.    {
  187.       for (i=0;i<n;i++)
  188.       {
  189.      x[i] /=n;
  190.      y[i] /=n;
  191.       }
  192.    }
  193.  
  194.    for(i=0;i<(n-1);i++)
  195.    {
  196.       if(i<j)
  197.       {
  198.      temp=x[i];
  199.      x[i]=x[j];
  200.      x[j]=temp;
  201.      temp=y[i];
  202.      y[i]=y[j];
  203.      y[j]=temp;
  204.       }
  205.       k=n/2;
  206.       while(k<=j)
  207.       {
  208.      j -=k;
  209.      k /=2;
  210.       }
  211.       j +=k;
  212.    }
  213.    maxpower=0;
  214.    i=n;
  215.    while(i>1)
  216.    {
  217.       maxpower++;
  218.       i /=2;
  219.    }
  220.    harm=2*PI/n;
  221.    for(i=0;i<n;i++)
  222.    {
  223.       sin_coef[i]=flag*sin(harm*i);
  224.       cos_coef[i]=cos(harm*i);
  225.    }
  226.    a=2;
  227.    b=1;
  228.  
  229.    for(j=0;j<maxpower;j++)
  230.    {
  231.       pnt0=n/a;
  232.       pnt1=0;
  233.       for(k=0;k<b;k++)
  234.       {
  235.      i=k;
  236.      while(i<n)
  237.      {
  238.         arg=i+b;
  239.         if(k==0)
  240.         {
  241.            prodreal=x[arg];
  242.            prodimag=y[arg];
  243.         }
  244.         else
  245.         {
  246.            prodreal=x[arg]*cos_coef[pnt1]-y[arg]*sin_coef[pnt1];
  247.            prodimag=x[arg]*sin_coef[pnt1]+y[arg]*cos_coef[pnt1];
  248.         }
  249.         x[arg]=x[i]-prodreal;
  250.         y[arg]=y[i]-prodimag;
  251.         x[i] +=prodreal;
  252.         y[i] +=prodimag;
  253.         i +=a;
  254.      }
  255.      pnt1 +=pnt0;
  256.       }
  257.       a *=2;
  258.       b *=2;
  259.    }
  260.    free(cos_coef);
  261.    free(sin_coef);
  262. }
  263.  
  264.  
  265.  
  266.  
  267.  
  268. /* get_filter() builds a high_pass or low_pass filter.  The arguments are: */
  269. /* FILTER is the array for the filter coefficients (frequency-domain */
  270. /* AMPLITUDE coefficients [0 to 1]), N is the number of elements in FILTER */
  271. /* SIGN is -1 for lowpass, +1 for highpass; NYQUIST is the nyquist freq */
  272. /* for the FILTER data set, CUTOFF is the filter cutoff frequency */
  273. /* ROLLOFF is the rolloff in dBells per filter_unit; */
  274. /* FILTER_UNIT is 2.0 for rolloff as dB/octave or 10.0 for dB/decade; */
  275. /* get_filter() will create positive/negative format amplitude filter */
  276. /* coefficients ONLY. The function returns the coefficients in FILTER */
  277.  
  278. double * get_filter(filter, n, sign, nyquist, cutoff, rolloff, filter_units)
  279. double filter[];
  280. int n, sign;
  281. double nyquist, cutoff, rolloff, filter_units;
  282. {
  283.     double freq_interval, active_freq;
  284.     double half_power, exp_var, octaves;
  285.     double filter_factor;
  286.     int i,m;
  287.  
  288.     half_power=sqrt(0.5);
  289.     filter_factor = log10(filter_units);
  290.     m=n/2;
  291.     freq_interval=nyquist/m;
  292.  
  293.     if(sign== (-1))
  294.     filter[0]=1.0;
  295.     else
  296.     filter[0]=0.0;
  297.  
  298.     for(i=1;i<=m;i++)
  299.     {
  300.     active_freq=freq_interval*i;
  301.  
  302.         /* number of octaves: */
  303.     octaves=log10(active_freq/cutoff)/filter_factor;
  304.         /* 20=constant for dB formula as amplitude */
  305.     exp_var=rolloff*octaves*sign/20.0;
  306.     exp_var=pow(10.0,exp_var);
  307.     filter[i]=exp_var*half_power;
  308.     if(filter[i]>1.0)
  309.         filter[i]=1.0;
  310.     if(filter[i]<MINVAL)
  311.     filter[i]=0.0;
  312.     if(sign==(-1) && filter[i]==0.0)
  313.         for( ;i<=m;i++)
  314.         filter[i]=0.0;
  315.     if(sign==1 && filter[i]==1.0)
  316.         for( ;i<=m;i++)
  317.         filter[i]=1.0;
  318.     }
  319.     for(i=1;i<m;i++)
  320.     filter[m+i]=filter[m-i];
  321.     unscramble_transform(filter,n);
  322.  
  323.     return filter;
  324.  
  325. }
  326.  
  327.  
  328.  
  329.  
  330.  
  331. /* int_scale() is a routine to scale integer values to fit graph */
  332. /* moved from herc.c 10/89 by T Clune */
  333.  
  334. int int_scale(in_val, in_min, in_max, out_min, out_max)
  335.  
  336.         /* NOTICE THAT ALL PARAMETERS ARE INT */
  337.         /* in_val = number to be scaled */
  338.         /* in_min = smallest possible in_val */
  339.         /* in_max = largest possible in_val */
  340.         /* out_min =